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High order ring-exchange interactions are crucial for the study of quantum fluctuations on many 
highly frustrated systems. A versatile and efficient quantum Monte Carlo method, which can handle 
finite and essentially zero temperature and canonical and grand-canonical ensembles, has long been 
sought. In this paper, we present the first exact quantum Monte Carlo study of a model of hard-core 
bosons with sixth order ring-exchange interactions on a two-dimensional kagome lattice. By using 
the Stochastic Green Function algorithm with global space-time update, we show that the system 
becomes unstable in the limit of large ring-exchange interactions. It undergoes a phase separation 
at all fillings, except at 5 and § fillings for which the superfluid density vanishes and an unusual 
mixed valence bond and charge density ordered solid is formed. This explains the universal features 
seen in previous studies on various different models, such as the transverse field Ising models, on a 
kagome lattice near the classical limit. 
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Interest in ring-exchange interactions in quantum 
many-body systems has a long history originating from 
the study of quantum solids, a typical example being 
solid Helium-3 [l[ 0. Recently the study of ring- 
exchange interactions has resurged with boson and spin 
models In particular, multiple particle exchange 

has been suggested as a possible candidate to induce a 
normal "Bose metal" or "Bosc liquid" phase at zero tem- 
perature, in which there are no broken symmetries associ- 
ated with superfluidity or charge density wave phases [§- 
Studies on a square lattice with four-site ring- 
exchange plaquettes suggest that ordered phases always 
dominate |4|-|2j. However, these results still leave the pos- 
sibility that a Bosc liquid or spin liquid phase may exist 
in certain frustrated lattices 



aii2|. 



Here we explore a frustrated kagome lattice built from 
corner sharing triangles in two dimensions. The most 
important characteristic of highly frustrated quantum 
systems is their ground state degeneracy in the classi- 
cal limit 13 , [3] ■ In such systems the classical ground 
state manifold, which is given by the least frustrated Ising 
spin configurations, results in local constraints in each 
unit. One of the prominent problems in frustrated quan- 
tum magnetism is whether quantum fluctuations par- 
tially alleviate the classical ground state degeneracy via 
the order-by-disorder mechanism or by forming aquan- 
tum spin liquid driven by quantum fluctuations [14|-|l6[. 



A natural route to study this problem is to construct an 
effective theory for the low energy quantum fluctuations 
which is confined to the Hilbert space of the degenerate 
classical ground state manifold. From the strong cou- 
pling perturbation theory, the leading term for the quan- 
tum fluctuations involves multiple-spin loop flips which 
form the ring-exchange term. For kagome and pyrochlore 
lattices, the smallest ring-exchange term appears at the 
sixth order [i| [ll|, [13] ■ Therefore, the study of quantum 
fluctuations on some of the most important highly frus- 



trated lattices naturally involves effective models with 
multiple-spin ring-exchange. 

Recently those ring-exchange models have been formu- 
lated in terms of gauge theories 0, [T3, EH ■ The gauge 
theory formulation is essentially the manifestation of the 
constraint for the projected Hilbert space within the clas- 
sical degenerate manifold, and the effective gauge theory 
is, in turn, studied by a duality mapping. Various phases 
have been proposed for different models based on this 
type of calculations. 

It is extremely important to have a systematic un- 
biased numerical method to test the various proposals 
for exotic ordered valence bond and disordered spin liq- 
uid phases. However, numerical studies on those models 
have proved to be rather difficult. There are some re- 
cent studies using the World-line and the Stochastic Se- 
ries Expansion algorithms with loop update schemes on 
four-site ring-exchange models 0-01 • F° r these models, 
ring-exchange interactions are described by a term that 
performs a correlated hopping of two particles in opposite 
directions, with no contribution to the winding. Thus 
one can expect the superfluidity to be destroyed when 
this term is dominant. However, instead of the proposed 
Bosc liquid phase, these exact studies have found charge 
density waves, valence bond solids, or phase separation. 

The model we consider in this work consists of hard- 
core bosons on a two-dimensional kagome lattice (Fig.[l}. 
The underlying Bravais lattice is spanned by the basis 
vectors (ai,a2) with lengths chosen as unity, and the 
kagome lattice is obtained by duplicating a set of three 
sites Si(|,0), S 2 {0,1), S 3 (±, |) (gray symbols). The re- 
ciprocal lattice is spanned by the vectors (61,62), each 
with length ^Txj\fZ. 

The Hamiltonian takes the form (we use periodic 
boundary conditions) 



H 



H.c. 



2 




FIG. 1: (Color online) The kagome lattice and the effect of 
the different terms in the Hamiltonian. Each site is shared by 
two hexagons, so the total number of sites is three times the 
number of hexagons (for a periodic lattice). The usual kinetic 
term t allows the particles to hop between near-neighboring 
sites. The ring-exchange term K performs a correlated hop- 
ping of three particles within the same hexagon. This process 
is possible only if the hexagon contains exactly three non- 
near-neighboring particles. The figure also shows our conven- 
tion for the labels of the sites of a given hexagon. 
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where a] and a i are the creation and annihilation op- 
erators of a hard-core boson on site i. These opera- 
tors satisfy fermionic anti-commutation rules when act- 
ing on the same site, of = 0, a\ 2 = 0, {a if o|} = 1, 
and bosonic commutation rules when acting on differ- 
ent sites, [a^aj = 0, [aj,aj] = 0, [a 4 ,a]] = 0. The 
usual kinetic term allows the particles to hop between 
near- neighboring sites i and j. A sextic ring-exchange 
term allows three particles to perform a correlated hop- 
ping within the same hexagon, with the over all 
hexagons. We study the model as a function of filling 
factor and K/t. 

This model can be mapped to the U(\) lattice gauge 
theory with a softened hard-core constraint, which plays 
a crucial role in the study of quantum fluctuations in 
some highly frustrated models 0, [ll[ . Notably, similar 
models can possibly be realized in optical lattices with 
dipolar bosons (l9T |. Since the sextic ring-exchange term 
couples six different sites at a time, an analytical treat- 
ment of the Hamiltonian (TTJ) is rather complicated 
On the other hand, most current quantum Monte Carlo 
(QMC) methods used on a square lattice 0-0] need to 
decompose the Hamiltonian as a sum of two-site coupling 
terms, which is not possible in our case. These methods 
need such a decomposition because they use a loop up- 
date that consists of building a closed path in space and 
imaginary time, and raising or lowering the occupation 
number of the sites that are visited. By construction, 
these loops can update only one or two sites at a time, 
and are therefore suitable only for one-site or two-site 
coupling terms. In brief, the treatment of fourth order 
ring-exchange terms in previous studies required some 
special developments 0, 0] that are not easy to general- 



ize to sixth or higher order. 

However, the Stochastic Green Function (SGF) 



rithm [2fj, |2l| with global space-time update 



not make use of such loop updates, instead, it performs 
a direct sampling of the partition function by distribut- 
ing Hamiltonian terms randomly in space and imaginary 
time. No particular decomposition of the Hamiltonian is 
required. Therefore the update procedure is totally in- 
dependent of the structure of the Hamiltonian, and can 
be applied to any n-site coupling term with efficiency. In 
this work, we use the SGF algorithm and perform simu- 
lations of systems with sizes up to 18 x 18 hexagons (972 
sites). 

While the SGF algorithm is designed to work in the 
canonical ensemble (CE), an extension [22j allows us to 
simulate the grand-canonical ensemble (GCE). In the fol- 
lowing we take advantage of this flexibility. We find con- 
venient to use the GCE for studying the stability of the 
system. However, it is much easier to use the CE on pa- 
rameter regions where the system is stable only at spe- 
cific densities. This is, to the best of our knowledge, the 
first time that an unbiased QMC method allows the sim- 
ulation of sixth order coupling terms in both the CE and 
GCE at both finite and essentially zero temperature. For 
the GCE we add the usual term — p,J\f to the Hamilto- 
nian ([T]) , where p, is the chemical potential and N is the 
total particle number operator. In the CE \x is not a con- 
trol parameter, but it is measured at zero temperature 
as (i(N) = E(N + 1) - E(N), where E = (H) and N is 
the controlled number of particles. 

In the following, we show that the model contains va- 
lence bond solids at densities p — | and p= \ when the 
ring-exchange term is dominant. For our study, we will 
consider the supcrfluid density 



) 



up 



(2) 



where (3 is the inverse temperature and W\ and W% 
are the winding numbers measured in the two direc- 
tions ti\ and tii- We also consider the structure factor 
S(k) = (|n(fc)| 2 ), with 

1 
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where L is the linear system size (counting hexagons), p 
is the density, h p = a p a p is the number operator on site 
p, and f p is the position of site p. Note that the subtrac- 
tion of p in ((3]) is meant to get rid of the lattice Bragg 
peaks. In order to study the ground state properties, we 
systematically use the inverse temperature ft = 2L/t. 

We start by illustrating in Fig. [2] the exactness of the 
SGF results, by comparing them with an exact diagonal- 
ization on a lattice with 2x2 hexagons (12 sites) and six 
particles. The figure shows an excellent agreement for 
the energy E and the supcrfluid density p s . 
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FIG. 2: (Color online) Comparison between an exact diago- 
nalization and the SGF algorithm for the energy E and the 
superfluid density p„. 



It is useful to analyze first the stability of the system. 
This is easily done by looking at the total density p as 
a function of the chemical potential p. The slope of this 
function is proportional to the compressibility of the sys- 
tem, so an instability of the system results either in a 
negative slope in the CE, or a discontinuity of the total 
density in the GCE. We first discuss the small K region, 
and search for any phase transition as we increase the 
value of K. Figure [3] shows results for K = 5t and a 
system with 4x4 hexagons (48 sites) in the CE and 
12 x 12 hexagons (432 sites) in the GCE. The agreement 
between the two curves reveals that finite size effects are 
sufficiently small to ensure the equivalence of the two 
ensembles. For this strength of the ring-exchange inter- 
actions the slope of the curve remains positive and finite 
at all fillings, so there is no instability. 
1 



— 4x4 hexagons (48 sites) - Canonical ensemble 
• 12x12 hexagons (432 sites) - Grand-canonical ensemble 




FIG. 3: (Color online) The total density p as a function of 
the chemical potential p for K = 5t. Results for a lattice 
with 4x4 hexagons in the canonical ensemble and 12 x 12 
hexagons in the grand-canonical ensemble are shown. The 
agreement between the two curves reveals that finite size ef- 
fects are sufficiently small to ensure the equivalence of the two 
ensembles. 

The situation changes as K increases. Fig. [4] shows 
the total energy E as a function of the density p in the 
CE for K/t = 10 and K/t = 25. We show only data 
for p < since the data for p > \ can be deducted by 
particle-hole symmetry. The cyan dotted line serves as 
a guide to indicate that the curvature of the K/t = 10 



curve is positive for p £ [|; |1 , while the orange dotted 
lines emphasize regions with negative curvatures. Sys- 
tems for which the energy has a negative second deriva- 
tive with respect to the density are thcrmodynamically 
unstable and undergo a phase separation. We also note 
the presence of a kink at p = i for both curves, which in- 
dicates a gapped phase. We conclude that for K/t = 10 
the system is compressible for p £ ] i ; | [ and in a solid 
phase for p = i and p = |, while it is unstable for all 
other densities. For K/t = 25 the system is unstable 
for all densities, except for p = -| and p = | for which 
the phase is solid. The insets correspond to GCE results 
that are in agreement and which, in addition to showing 
"jumps" over the "forbidden" densities, show the exis- 
tence of a "numerical hysteresis" , depending on if the 
initial state is empty (ascendant hysteresis) or full (de- 
scendant hysteresis). This phenomenon is observed in 
many Monte Carlo studies 0, [24j, including the four- 
site ring-exchange model Q. For these parameters the 
grand-canonical algorithm is no longer able to sample all 
contributing states, which is an indication that the sys- 
tem undergoes a spontaneous symmetry breaking in the 
thermodynamic limit. This also supports the importance 
of the ability of the SGF algorithm to work in the canon- 
ical ensemble. 
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FIG. 4: (Color online) The energy as a function of the density 
for K/t = 10 and K/t = 25 in the CE. The insets show for 
comparison the density as a function of the chemical potential 
in the GCE, with a numerical hysteresis. See text for details. 



Since only the fillings P = \ and p = | are stable in 
the large K limit, it is convenient to analyze how the su- 
perfluid density p s is destroyed as a function of K/t by 
working in the CE. We show in Fig. [5] results for p = \ 
only, the case with p = | being identical because of the 
particle-hole symmetry. There exists a critical value of 
the ring-exchange interaction, K c sa 9t where the super- 
fluid density vanishes. We show results for 6 x 6 hexagons 
(108 sites) and 12 x 12 hexagons (432 sites) to illustrate 
that finite size effects are small. 

For t K, the model is in the free hard-core limit, 
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FIG. 5: (Color online) The superfluid density p a as a func- 
tion of K/t for p — 1. The superfluid density is completely 
destroyed at the critical value K c « 9i. 



thus the superfluid phase prevails against solid or va- 
lence bond ordering. In the opposite t <C K limit, we 
expect the ground state to be formed by configurations 
that maximize the number of resonating hexagons. Be- 
cause a resonance can occur only if the hexagon contains 
exactly three non-first-neighboring particles, it is not pos- 
sible to have two resonating neighboring hexagons at first 
order in K. This suggests that a phase with a maximum 
number of resonating hexagons can be formed at p = | 
by having all resonating hexagons surrounded by empty 
hexagons (Fig. [S] left), or at p = | by having all res- 
onating hexagons surrounded by hexagons that contain 
three localized particles on the vertices of the triangles 
that are not shared with the resonating hexagons. (Fig. [6] 
right). These phases can be understood as a form of va- 
lence bond solid in which each valence bond now involves 
six sites, in contrast with the usual two-site singlets. We 
note that the phase has been suggested in Kagome lat- 
tice on other models, such as the transverse field Ising 
model [HHH. 




FIG. 6: (Color online) Two possible ordering patterns for 
p = \ (left) and p = | (right). The circles inside the hexagons 
represent three resonating particles. The red dots represent 
localized particles. We note that both orderings have three- 
fold degeneracy. These two uniform ordering patterns are 
only formed for p = | and p = |. Other fillings are unstable 
towards phase separation. 



correspond to the minima of the intensity (black and 
blue). For K = 5i (middle panel), high intensity regions 
(red and yellow) appear in the center of triangles formed 
by the Bragg peaks. In the insulating phase, K = 15t 
(right panel) , the high intensity regions become localized 
at k max = (|, i) in the (61,62) basis (and symmetry re- 
lated momenta), and form a honeycomb lattice. Since 
k m ax is parallel to a,\ and ||fc TOaa; || = 4^, the insulating 
phase has features that appear with a spatial frequency 
of I in the a\ direction, in agreement with Fig. [6] 




FIG. 7: (Color online) The structure factor S{k) for 18 x 18 
hexagons (972 sites), p = | and K = (left panel), K = 5t 
(middle panel), and K = 15t (right panel). The black re- 
gions correspond to the locations of the lattice Bragg peaks. 
As K increases, peaks with maximum intensity develop at 
k = + ^62 and symmetry related momenta. 



To conclude, we study hard-core bosons on a kagome 
lattice with a sextic ring-exchange term using the 
Stochastic Green Function algorithm Wc find 

that the system becomes unstable towards phase sepa- 
ration as the ring-exchange interaction increases, except 
for densities p — \ and p — | where the superfluid is de- 
stroyed. Here we observe an unusual valence bond solid 
phase with resonances involving six sites simultaneously. 
The hysteresis obtained in the quantum Monte Carlo in- 
dicates that the phase transition between the superfluid 
and the valence bond solid at p = | and p = | is first 
order. Models with higher order ring-exchange played an 
important role in the search of spin liquid phases due to 
their relation with the gauge theory deconfined phase. 
This work showcases the power of the Stochastic Green 
Function algorithm to study models with higher order 
ring-exchange interactions which have hitherto been very 
challenging for other Monte Carlo methods. 
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This scenario is confirmed by looking at the structure 
factor. Figure [7] shows S(k) for 18 x 18 hexagons (972 
sites) with p = §, K = 0, K = 5t, and K = 15*. For 
K = (left panel), only lattice Bragg peaks arise and 
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